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Abstract. We develop a systematic metiiod to obtain the solution of the colUsionless Boltzmann equation which 
describes the growth of large-scale structures as a perturbative series over the initial density perturbations. We 
give an explicit calculation of the second-order terms which are shown to agree with the results obtained from the 
hydrodynamical description of the system. Then, we explain that this identity extends to all orders of perturbation 
theory and that the perturbative series actually diverge for hierarchical scenarios. However, since the coUisionless 
Boltzmann equation provides the exact description of the dynamics (including the non-linear regime) these results 
may serve as a basis for a study of the non-linear regime. In particular, we derive a non-perturbative quadratic 
integral equation which explicitly relates the actual non-linear distribution function to the initial conditions 
(more precisely, to the linear growing mode). This allows us to write an explicit path- integral expression for the 
probability distribution of the exact non-linear density field. 
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1. Introduction 

In usual cosmological models large-scale structures in the 
, universe have formed through the growth of small ini- 
' tial density perturbations by gravitational instability, see 
pccbles (1980)| . Moreover, in most cases of cosmological 
■ interest the power increases at small scales as in the CDM 
, model (Peebles (1982)). This leads to a hierarchical sce- 
nario of structure formation where smaller scales become 
non- linear first. Once they enter the non-linear regime 
high-density fluctuations form massive objects which will 
give birth to galaxies or clusters. Thus, it is important to 
understand the dynamics of these density perturbations 
in order to describe the large-scale structures we observe 
in the present universe. 

The distribution of matter is usually described as a 
pressure-less fluid through the standard equations of hy- 
drodynamics (continuity and Euler eq.) coupled to the 
Poisson equation for the gravitational potential. Then, one 
often looks for a solution of the equations of motion un- 



in the non-linear regime shell-crossing should play a key 
role (e.g., through virialization processes). 

Therefore, in order to obtain a complete de- 
scription of the dynamics one needs to go beyond 
the standard approximation of a pressure-less fluid. 
Then, one might add a pressure term to the hy- 
drodynamical equations in order to describe the ef- 
fects of multi-streaming, which provide some ve- 
locity dispersion (e.g., Buchert & Dominguez (1998) 



Buchert et al. (1999)). However, this is not very satisfac- 



der the form of a perturbative expansion (e.g.. Fry (1984) 
Goroff ct al. (1986), Bernardeau (1994)|). In fact, one can- 



not hope to describe the non-linear regime within such a 
framework since this hydrodynamical description becomes 
inexact as soon as shell-crossing appears. Then, there is 
no longer only one velocity at each point. Note that for 
hierarchical scenarios with no small-scale cutoff (i.e. the 
amplitude of density fluctuations diverges at small scales 
i? ^ 0) shell-crossing occurs as soon as i > 0. Moreover, 



tory because coUisionless and gaseous systems can ex- 
hibit very different behaviours starting from the same 
initial conditions. For instance, some physical processes 
are specific to coUisionless systems (e.g., Landau damp- 
ing, phase-mixing) and stability conditions for similar col- 
lisionless and gaseous systems are not always the same, 
see Binney fc Tremainc (1987) , Hence it is important to 
investigate the dynamics of density fluctuations within the 
framework of the coUisionless Boltzmann equation, which 
provides an exact description of the system both in the 
linear and non-linear regimes. 

Since this is a rather difficult task we may first try to 
devise a perturbative approach. This may then serve as 
a basis for further developments to tackle the non-linear 
regime. Thus, in this article we build a perturbative ap- 
proach to the coUisionless Boltzmann equation. One could 
also start from the BBGKY hierarchy, which can be de- 
rived from the Boltzmann equation. Such a study was per- 
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formed in Bharadwaj (1996) where it was shown that the 
lowest non-linear order matches the results of the usual hy- 
drodynamical approach. We recover this result here and 
we explain that indeed a perturbative method cannot han- 
dle shell-crossing. However, in the building of perturbative 
theory we also derive in this article some non-perturbative 
results which may prove to be useful for a study of the 
non-linear regime. 

This paper is organized as follows. In Sect.y we write 
the equations of motion which describe the dynamics of 
the collisionless fluid. Then, in Sect.y we develop a sys- 
tematic procedure to obtain the fluctuations of the dis- 
tribution function as a perturbative series over the initial 
perturbations (or more precisely over the linear growing 
mode) in the case of a critical density universe. In doing so 
we also derive a non-perturbative quadratic integral equa- 
tion which directly relates the actual non-linear distribu- 
tion function to this linear mode. Next, we give in Sect.H 
the explicit calculations of the second-order terms. We ex- 
plain in Sect.H why we recover the results of the hydrody- 
namical approach and why the perturbative series must 
diverge for hierarchical scenarios (with no small-scale cut- 
off). Then, we show in Sect.O how to extend this method 
to arbitrary cosmologies. Finally, we derive in Sect.H an 
explicit path-integral expression which describes the sta- 
tistical properties of the exact non-linear density field. 



2. Equations of motion 

The gravitational dynamics of a fluid of collisionless par- 
ticles of mass m is given by the collisionless Boltzmann 
equation ( Peebles (1980) ), coupled with the Poisson equa- 
tion: 



( df ^ p df ^d^ df _ ^ 

dt ma? dx 9x dp 

AttQ , _, 

A0= [P-Pj , p[x,t}^m 

a 



/ /(x,p,t) dp 



(1) 



where /(x, p, t) is the distribution function (phase-space 
density), (j) is the gravitational potential, p(x, t) is the co- 
moving density, p is the mean comoving density of the uni- 
verse (hence it is independent of time), a{t) is the scale 
factor and the momentum p is related to the comoving 
coordinate x of a particle by p = ma?±. Thus, eq.(P 
describes the growth of gravitational structures in an ex- 
panding universe which is homogeneous on large scales. 

Note that in order to derive the collisionless Boltzmann 
equation one needs to neglect the two-particle correla- 
tions. Then, the gravitational interaction is described by 
a smooth gravitational potential which neglects the ef- 
fects of the discrete character of the distribution of mat- 



ter (e.g., Peebles (1980) ). This becomes exact in the con- 
tinuous limit TO — > which we perform below in eq.(y). 
However, one must also add some small-scale cutoff to 
the power-spectrum of the density fluctuations in order to 
guarantee small-scale smoothness. 



The definition itself of the problem we investigate here 
also involves an implicit "regularization" of the gravita- 
tional interaction. Indeed, within an infinite uniform uni- 
verse the gravitational force is not well defined since the 
pull due to the matter located within a small solid angle 
dfi diverges as we integrate over all shells up to infinity. 
This problem is cured by first integrating the relevant in- 
tegrals over angles, and next over radial distances (e.g., 
Peebles (1980)). This can also be seen as the introduction 
of a large-scale cutoff L for the gravitational interaction 
which preserves rotational symmetry. More precisely, this 
means that in Fourier space we integrate on the wavenum- 
ber k over |k| > k^ and in the final results we take the 
limit fee ^ 0. In fact, except in eq.(|7q) in Sectj^ this fea- 
ture will not show up in the equations we encounter in 
this article so that we can directly take kc = 0. 

Next, we can absorb the mass m into the distribution 
function / and the momentum p: / ^ //tti'*, p —>■ rnp, 
and we obtain: 



a/ ^ P_ 9/ _ 90 a/ ^ Q 

dt a^ dx 9x dp 



A<j>-- 



AttG 



{P~P) , P(x,t) 



y/(x,p,t)dp 



(2) 



which is also valid in the continuous limit ?7i — > 0. Here 
/(x, p, t)d^x(l^p is the mass enclosed in the phase-space 
element A^xA^p and the momentum p verifies: 



a X. 



(3) 



In order to single out the deviations of the distribution 
function from the case of a perfectly homogeneous universe 
we define the perturbation 6f(x, p, t) by: 



/(x, p,t)=p Snip) + p (5/(x, p, t) 



(4) 



where So is Dirac's function (here in 3 dimensions). Then, 
the density contrast (5(x, t) is simply given by: 



S{x,t)^ J Sf{x,p,t) dp. 



(5) 



From the definition 1^ and the equation of motion (g) we 
obtain the evolution equation for Sf: 



dt a^ dx dx dp dx 9p 

A(j)^^pS{x,t) 
a 



(6) 



Note the term ^^(p) which is sometimes forgotten in the 
literature. Next, we define the spatial Fourier transform 
of the field 6f by: 



^/(x,p,t) = ydke*-^J/(k,p,t) 

5/(k, p, t) = -i^ J dx e-''^'^ 6 fix, p, t) 



(7) 
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To simplify the notations, we note the distribution func- 
tion by Sf, both in real space x and in Fourier space 
k. Its argument removes any ambiguity. We also define 
the Fourier transforms (5(k) and 0(k) of the density con- 
trast and of the gravitational potential as in eq. {uh . Then, 
Poisson's equation reads: 



/fc20(k) = —p (5(k) 



(8) 



After substitution of eq.(|g) into the equation of motion 
(^ we obtain: 






dt 



-f.l^/dk'dp'^/(k',p')^.^(k-k',p)=0. (9) 

This non-linear and non-local equation describes the dy- 
namics of the coUisionless fluid. 



3. Building a perturbative expansion 

3.1. General framework 

Here we consider the case of a critical density universe 
r^m — 1. Then, the average comoving density is given by: 



Qirgt^ 



(10) 



We normalize the time scale such that the Hubble time 
today io (i-e. at z = 0) satisfies io = Ij while the scale 
factor is a = 1. Then, we define a new time coordinate t 
by: 

r=iln(tAo)-^lnt, t = e'^ a = e"^ (11) 

and eq.(||) reads: 

^ + 3z(k.p)e-^5/ + 2z e^ii.^(p) j dp' ^/(k, p') 
+2z e- j dk'dp'5/(k', p')^.^(k - k', p) = 0. (12) 

This equation is quadratic over the distribution function 
5f. In order to make apparent the structure of this equa- 
tion of motion it is convenient to separate the linear and 
quadratic terms. Thus, we can also write eq.(|l2|) as: 



and: 



iir(r;ri,r2) = 2i e^ Soiji - t)6d{t2 - t) 

X(5D(ki + k2 - k) p-.-^(p2 - p). 



(15) 



In the r.h.s. in eq.(13) the integration over ri and Vi runs 
over all space. In particular, the time coordinates t\ and 
Ti are integrated over the whole range ] — oo, -|-(X)[. Note 
however that the Dirac factor 5y)(j\ — t)5d{t2 — t) which 
appears in the kernel K implies that the r.h.s. only de- 
pends on the distribution function 5f at the time r, as in 
eq.(|l^). This also ensures that causality is not violated. 
In the relation (|l3| ) we also introduced a "coupling con- 
stant" g which is equal to unity: g = \. The advantage 
of generalizing the quadratic equation (13) to arbitrary g 
will appear later on. 

3.2. Perturbative expansion 

The dynamics of the distribution function 5f is fully de- 
termined by the non-linear equation (|l3| ) supplemented 
by initial conditions. At early times r — > — oo the uni- 
verse becomes homogeneous and the velocity is given by 
the Hubble flow (i.e. p ^ 0) so that the perturbations 
vanish: Sf ^ 0. Thus, at early times (or at large scales 
for CDM-like initial conditions) the equations of motion 
can be linearized and the distribution function is equal to 
the linear solution 77(r). As time goes on, the perturba- 
tion grows and it eventually enters the non-linear regime. 
In the mildly non-linear regime, one can seek a solution of 
eq.dl3) in the form of a perturbative expansion over the 
linear mode ri{r). Thus, we write the expansion: 



<5/(r) = £<5/(")(r) 



with: 



5/(1) (r) EE 77(r) and (5/(")(r) - r^''' 



(16) 



(17) 



That is, the term of order n is of the order of rj to the 
power n. Substituting the expansion (16) into the equation 
of motion (|l^) we obtain: 



0.77 = 

and for n > 2: 



(18) 



/• , , 0.(5/(")=5 /dridrzi^ y (5/(")(ri),5/("-™)(r2). (19) 

{0.5f){r)=gjd'r,d''r2K{r;r^,r2)Sfir,)Sfir2) (13) 'J ^, 



where we note r the 7-dimensional coordinate r — (k, p, r) 
and we introduced the linear operator O and the kernel 
K defined by: 



(0.5/)(r)^^+3z(k.p)e--5/ 



(14) 



Thus, the linear mode /^(r) must belong to the kernel of 
the operator O. Since we require 77 ^ this implies that 
O is not invertible. We shall check below that indeed the 
kernel of the operator O defined by eq. (M[) is not reduced 
to {0}. Then, the higher-order terms (5/*^") are given by 
the recursion ([l9|). 

However, in order to obtain an explicit expression for 
(5/("-) we must "invert" the operator O in eq.(|ig|). More 
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precisely, we need to build a kernel K{r; ri, r2) which sat- 
isfies: 



(O.AT(r;ri,r2) = X(r;ri,r2). 



(20) 



Indeed, if such a kernel K exists, we can solve the recursion 
(19) by choosing for the term of order n the value: 



5f 



in) 



9 



dndrsK ^<5/('")(ri)5/("-'")(r2). (21) 



This relation fully defines the expansion (|l^): it provides 
an explicit procedure to compute any order term 5f^'^\ 
in a recursive fashion. Note that it is not obvious a priori 
that a solution K to cq.(|2C|) exists. On the other hand, 
if a solution exists it is not unique since we can add to 
any peculiar solution K any factorized kernel of the form 
77(r)iJ(ri,r2) where H is arbitrary, as shown by eq.(p^. 
This amounts to add a factor -q to all terms Sf^^K We 
shall come back to this point below. 

Thus, we now need to construct a solution K to 



eq.(20). To do so, we use the following trick. First, we 
define a new operator O^i by: 



O-i(r) 



dae--°W. 



(22) 



Since the operator 0_i is only an intermediate tool eq. (p2[) 
should be interpreted in a formal sense. In particular, one 
may introduce a large cutoff A in the intermediate steps 
of the calculation for the integration over a in the defini- 
tion (p3) and take the limit A ^ oo in the final results. 
The reason behind the definition (E2) is that this operator 
O-i is formally the "inverse" of the operator O. Indeed, 
expanding the exponential e~'^ we can write: 



(^"'^''■/)=E /„!), 0'^..f = -0.e-^°.f (23) 



da 



where /(r) is an arbitrary function of r. Using eg. (|23|) we 
obtain: 



O.O-i.f ^ I daO-e-"^./ 



(24) 



if all integrals converge. Here, we must point out that 
eq.(24) is not a rigorous proof. In particular, using the 
same arguments we could also obtain O^i.O.f = f. 
However, we know that O^i.O cannot be the identity op- 
erator since the kernel of O is not reduced to {0}, as shown 
by eq.(nS) which implies O-i.O.rj = 0. Nevertheless, we 
can still work with the operator O-i. In particular, we 
define a kernel K{r; ri, r2) by: 



X(r;ri,r2)EEO_i(r).if(r;ri,r2). 



(25) 



If the operator O were invertible and all integrals were 
convergent the kernel K would be the unique solution of 
eq.(20). Here, the operator O is not invertible hence it is 
not obvious a priori that the kernel K defined by eq.(25) is 



a solution of cq.(EOh. However, if this procedure provides a 
finite kernel K we could expect that it obeys eq. (pO[) . We 
compute in AppM and App.H the kernel K defined by 
eq. (|25| ) . This yields the equivalent expressions (B.l) and 

we also explicitly check that 
is a solution to eq.(pO|). 



(B.3). As discussed in App. 

the kernel K written in eq.( p3.3D 

Then, we obtain the explicit perturbative expansion of 

the distribution function 6f through eq.(pl|). As noticed 

above, at each step n we could add to Sf^'^an element of 

the kernel of O like rj. However, since we require that (5/'"^ 

be of the order of rj to the power n (so that a perturbative 

expansion makes sense) this additional term is zero and 

the high-order terms 5/^"^ are indeed given by eq.(pi|) for 

n>2. 

Here we note that instead of the operator 0_i defined 
by eq.(^) we could have tried to use the operator 0_i 
defined by: 



/■oo 

d_i(r) = - / dae"^^'^ 



(26) 



Indeed, we can easily see that this operator also satisfies 
O.O-i = 1 in a formal sense, following the steps (|2^) and 
(P4|). However, only one of the two choices O-i and O-i 
is valid and leads to finite results while the other one in- 
volves divergent quantities. In fact, it is easy to see a priori 
that the correct choice is O-i from causality requirements. 
Indeed, from eq.(14) we note that O = d/dr + .... Then, 
since we have: 



q=0 ^' 



we can see that (^/^"•'(t) defined by eq.(21) using O-i 
involves the linear mode at earlier times r](T — a) (with 
(7 > 0), as required by causality. On the contrary, using 
O^i would imply that 6f^'^'>{T) depend on ?7(t + ct), that is 
on the linear mode at later times, which violates causality. 
As a consequence, the correct choice is the operator 0_i. 
However, the actual justification of cq.(p5|) is the explicit 
check (discussed in AppM) that the kernel K derived in 
Ap p.[A| and App.H from the definition (G3) is a solution to 
eq.(|2^. 

From the recursion ( |2l| ) we can also express (5/^"^ ex- 
plicitly in terms of the linear mode rj. Indeed, we can easily 
check that we can write eq.(pi|) as: 



Sf 



(») 



(r)-5"-i /"dri..dr„F„(r;ri 



77(ri)..77(r„)(28) 



for n > 1. Here we defined the kernels Fn by: 

i^i(r;ri) = (5z,(r-ri) (29) 

and for n > 2: 

i^„(r;ri,...,r„) = / dr'^r^ K{r;r[,r'^) 

n-l 

X ^ F,„(r'i;ri,..,rm)F„_„j(r2;r,„+i,..,r„) (30) 

m— 1 
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The expression (ESh clearly shows that the expansion 
is a perturbative expansion over powers of the linear mode 
77. It is also apparent from eq.(p8|) that rj should only con- 
sist of the linear growing mode. Indeed, the decaying mode 
would lead to divergences in eq.(|2q) as r ^ —00. This is 
not surprising since at early times the decaying mode be- 
comes increasingly large so that a perturbative approach 
fails. Of course, in principle it is possible to include de- 
caying modes if we set up the initial conditions at a finite 
time Ti (i.e. ti > 0). Note that in this case we cannot use 
the kernel K derived from eg. (p5|) . This can be seen from 
the fact that if we only include the growing mode but set 
up the initial conditions at i^ > we should recover the 
same results for t > ti. However, this cannot be achieved 
from eq. (|2l|) using the same kernel K since the integration 
range over ti and T2 is changed to ti > Ti and T2 > Ti. 

Fortunately, for practical purposes the decaying modes 
should play no role: ti ^ to so that the decaying modes 
become negligible before one enters the non-linear regime. 
Hence it is convenient to restrict 7y(r) to the linear growing 
mode and to take ti = 0. 

3.3. Non-perturbative integral equation 

Here, we note that the expression (gSJ) also shows that the 
expansion (|l^) can be formally reinterpreted as a pertur- 
bative expansion over powers of the "coupling constant" g. 
Thus, the recursive solution defined by eq.(17) and eq.(pi|) 
is also the perturbative solution of the integral equation: 



Sfi^) ='7(r)+5 



y dndrz i^(r;ri,r2) <5/(ri)<5/(r2) (31) 

where the perturbative parameter is g. Moreover, if we ap- 
ply to both sides of eq. (|3l|) the operator O, using eq.(po|), 
we obtain eq.(|l^). This means that the solution SJIt]] of 
eq. (|3l| ) is also a solution of the non-perturbative equation 
([l3|). Hence eq.(pl|) is actually equivalent to eq.([l3|) sup- 
plemented by the initial condition (5/ — > 77 at early times 
T -^ — cx). Indeed, eq.(13) admits only one solution (i.e., 
the distribution function Sf is completely specified by its 
value at some initial time ti since eq.(|l3|) involves a first- 
order time derivative). Thus, the problem of solving the 
collisionless Boltzmann equation with this initial condi- 
tion comes down to look for the solution Sf of eq. (|l]) for 
a given linear mode 77. We stress here that eq.(pl|) is non- 
perturbative, since it is equivalent to eq.dls). Note that in 
our case we shall eventually put g = 1. In particular, the 
relation ( pl| ) does not assume that the distribution func- 
tion Sf can be written as a perturbative series over the 
linear mode r], that is eq.(pih remains valid even if the 
perturbative series is only asymptotic (i.e. divergent). 

The eq.(^I]) is similar to the integral form of usual 
stochastic differential equations (e.g., Langevin equation) 
where rj would be called a noise term. Moreover, this for- 
mulation is very convenient since it reduces the three equa- 
tions we started from (Boltzmann eq., Poisson eq., initial 
conditions) to the only one eq. (|3l|) . Thus, the initial con- 
ditions r7(r) are encoded within the "evolution equation" 



itself, which fully determines the properties of the distribu- 
tion function Sf once the characteristics of 77(r) are spec- 
ified. Here we note that a quadratic equation of the form 
(pTl) was derived in ^coccimarro (2000) from the usual hy- 
drodynamical description. However, the latter formulation 
breaks down when shell-crossing occurs while eq.(|3l|) re- 
mains valid in the non-linear regime. 

3.4. Diagrams 

The perturbative expansion (E^) can be conveniently writ- 
ten as a sum of tree-diagrams with a three-leg vertex. To 
do so, we first define the symmetrized kernel Kg by: 



^s(r;ri,r2) 



- X(r;ri, 



r2) + i^(r;r2,ri) 



(32) 



Thus, the relations ( pi|) , ( pO| ) and (pl^) remain valid after 
we replace the kernel K by K^. Then, one can write the 
term (5/^") of order n as the sum over all trees built from 
(n — 1) three-leg vertices with n external points ri, ..,r„ 
over which we integrate. For instance, we show in FigJ^ 
the two diagrams which give the fourth-order term Sf^ '. 
To construct such diagrams one simply applies a recursive 
procedure of {n — 1) steps which builds the tree upward. 
To go from step k to step k + 1 one merely chooses an 
end-point (i.e. with no links upward) and connects to this 
point two "sons" through a three-leg vertex. At the first 
step one starts from the "root" r while the final end-points 
which are left are labeled ri, .., r„ (the order is irrelevant). 
Then, each vertex contributes a weight gKs while to each 
external point is affected a weight rjijCi). Finally, to obtain 
(5/(")(r) one simply integrates over the points ri, .., r„ and 
takes the sum over all diagrams. Note that each diagram of 
topology a needs to be multiplied by a multiplicity factor 



sf*"': 




(1) r 




Fig. 1. The two diagrams which give the fourth-order 
term Sf^^'{v). The filled circles correspond to the three- 
leg vertices, which contribute a weight gKg for each one, 
and the circles represent the external points ri, .., r„. The 
numbers (1) and (4) give the multiplicity a4_Q of each di- 
agram. 

It is instructive to consider the simple case where Sf 
and rj are mere numbers, rather than functions. Then we 
can absorb Kg into g and the integral equation ( pi| ) be- 
comes: 



Sf = V + 9 5f^ 



(33) 
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Let us note c„ the sum over all multiplicity factors associ- 
ated with the diagrams of order ji, that is: c„ = ^^ dn^a 
(in particular, C4 = 5). Then, the perturbative expansion 
M) reads: 



^/-E.' 



(34) 



with ci = 1. However, in this simple case we can directly 
solve the non- linear eq. ( p3| ) which yields (the choice among 
the two roots of the quadratic eq. ( pSJ ) is determined by the 
constraint Sf = r/ in the limit g ^ 0) : 



Sf 



1 - yi - 4gr? 
2ff 



(35) 



Then, expanding the square root which appears in (|35| 
we obtain the coefficients c„ : 



ci = 1 and for n > 2 



2(2n-3)! 
n!(n-2)!' 



(36) 



Note that these coefficients c„ are simply the Catalan 
numbers. From the expression (^5|) we can see that the 
perturbative expansion ( p4| ) diverges for \gr]\ > 1/4. Note 
that this occurs rather early as the second-order term is 
only one fourth of the first-order term. This means in 
particular that we cannot use a local approximation for 
eq.(31 ), which would consist of taking the mean of eq.(pl|) 
over a small phase-space volume V and approximating the 
integral in the r.h.s. by its value when ri and r2 are re- 
stricted to V. Of course, this result does not imply that 
the exact perturbative series diverges for large 77 (i.e. in 
the non- linear regime), however it suggests that such a 
behaviour would be quite likely. 

4. Explicit calculations up to the second-order 

We have described in the previous section the general 
structure of the evolution equation which governs the dy- 
namics of the distribution function Sf{r). In particular, 
we have shown how we can obtain a perturbative expan- 
sion of (5/(r) over powers of the linear growing mode ri{r). 
In this section, we apply these results to the specific case 
of gravitational clustering in an expanding universe (with 

n„, = 1). 



4.1. Linear growing mode 

First, we have to obtain the linear growing mode ri{r) (as 
explained in Sect. 3.2 we do not need the decaying modes). 
This function Tyfr) must belong to the kernel of the oper- 
ator O, see eq.([l8|), and grow with time. To get ?7(r) one 
can simply consider the linear growing mode of the hydro- 
dynamical approximation. Thus, in this approximation all 
particles at a given point x have the same peculiar velocity 
V defined by: 



In the linear regime the growing modes of the density and 
velocity fields are related by (Peebles (1980)): 



VL{k,T)=t^D+^6L{Kr)^td^SL{k,T) (38) 

where the growth factor D^{t) verifies D+{t) ex a{t) 
since we consider a critical density universe. The dot over 
"D+" or "d" stands for the time-derivative "dl?+/di" or 
"da/di". The subscript L refers to the "linear" regime. 
Thus, using eq.(|37|) we have in the linear regime: 



PL(k, t) ^iaa —^ (5L(k, r). 

Moreover, since ^^(k, r) ex a we can write: 

2 



5L(k,T) 



^2t 



(5Lo(k), PL(k, r) 



'r 



3r 



5w{k) 



fc2- 



(39) 



(40) 



Thus the linear growing mode which is solution of the 
hydrodynamical equations is: 

'5/hyL(^.P'^) = (l+'5L)(x,r)<5z,[p-pi(x,r)]-<5z5(p).(41) 

The Dirac distribution (5£)[p — Pl(x, t)] corresponds to 
the hydrodynamical description. However, one can check 
that this function 5fl^\^ is not a solution of eq.([18|). To 
obtain the correct linear mode ?7(r) which belongs to the 
kernel of O one simply needs to keep only the linear term 
in eq.([4l|). Thus, we write: 



(5/(^) (x, p, r) = 5l (x, t)5d (p) - Pl (x, t) 



85 



D 



In Fourier space we obtain for ry(r) 
?7(r) = 



9p 

5/(i)(r): 



P). (42) 



k 85 



D , 



<5Lo(k)fo(p) - ^-e'-5Loik)-.-^{p). (43) 



with r/ given by 



V = ax, hence p = av. 



(37) 



Then, one can easily check that O.rj 
eq.dl) andObyeq.(|l4l). 

Here, we can note that up to first-order our system can 
be described as a pressure-less fluid. Indeed, we have not 
included any velocity dispersion in the initial conditions. 
As is well-known, caustics will form as the density field 
evolves and multi-streaming regions appear. The pertur- 
bative approach should break down at this stage since it 
cannot go beyond these singularities. Thus, we shall check 
in the following sections that our perturbative treatment 
recovers the predictions of the standard perturbative ap- 
proach applied to the equations of motion derived within 
the hydrodynamic approximation (see Peebles (1980) or 
poroff et al. (1986)] for a presentation of this method). 
This clearly implies that the perturbative method can- 
not handle the multi-streaming regime. However, we stress 
that the equations of motion (1^), dl3) and ( pl| ) remain 
valid in the non-linear regime after multi-streaming ap- 
pears. 

In order to handle the multi-streaming regime in a per- 
turbative way, one may try to add a small velocity disper- 
sion to the initial conditions. Indeed, this will smooth out 
the caustics encountered for a pressure-less fluid. Such an 
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approach is investigated in Buchert & Dominguez (1998) 
and Buchert et al. (1999) for instance. However, it is clear 
that this is not sufficient to go beyond shell-crossing. 
Indeed, as non-linear structures form, shocks will ap- 
pear around virialized objects. This can be seen from 
the analytic solution of gaseous spherical collapse in 
Bertschinger (1985). Then, a perturbative approach can- 
not give any information on the inner region beyond the 
shock. In particular, one cannot rely on a simple equation 
of state of the form p — p{p) for the pressure. Moreover, it 
is well-known that coUisionless and gaseous systems can 
exhibit qualitatively different behaviours (e.g., existence 
of Landau damping and ph ase-mixing for the former, see 
Binney fc Tremaine (1987) ). This is why it is important 
to devise methods which are based on the coUisionless 
Boltzmann equation (|lj) . This article is a first step in this 
direction. 

4.2. Second-order term. Comparison with the 
hydrodynamical approach 

Once the linear growing mode 77 (r) is specified we can de- 
rive the higher-order terms from the recursion (|2|). Thus, 
using the expression (Q) and eq.( p3.l| ) obtained in AppJ^ 
for the kernel K we get after a simple calculation: 

(5/(2) (k, P,r)^ f dki <5Lo(ki)^Lo(k - ki) 



2k.(k-ki) 

3 |k-ki|2 



3 4.r kki 



7 



Snip) 



2i g^ki dSa f ^ 6i 5.rk-ki k dSo 
5^ l4'~d^ " 






kj k^ ■ dp ^^\ 



2 er^± 

9 kfdp 



k-ki 






To derive eq.(44) we used identities of the form: 



F(k.p) Sd{p) ^ FiO) Soip) 



(44) 



(45) 



and: 

k dSo , ■. 



_F'(0)^ Snip) 



(46) 



which are valid for any function F. This allows us to elim- 
inate the functions ip which appear in K in eq.( |B.l|) since 
we can express all factors in terms of ■0(s,O) — 1/s, see 
eq-dU). 

It is interesting to compare this result for 5f^'^' with 
the results obtained within the usual hydrodynamical ap- 
proach. In this latter case, the system is fully defined by 
the density contrast and the mean fiuid velocity at each 
point X (there is no velocity dispersion). Thus, we need to 
derive the density and the mean collective velocity from 
eq-dH). 



4.2.1. Density contrast 

From eq.(|5|) and eq.(^) we obtain the second-order term 
for the density contrast: 

(5(2) (k, r) = e*^ J dki 5Lo(ki)(5Lo(k - ki) 



'3k.ki 2k.kik.(k-ki) 

"" ' 7~kf^7 kj |k-ki|2 



(47) 



Going back to real space, we get after a change of vari- 
ables: 

,5(2) (x,t) = e^^y dkidk2 e'(ki+k2).x SLo{ki)6Lo{k2) 



ki.k2 

kf 



kik2 



(48) 



Thus, we recover the result obtained by the hydrodynam- 



ical approach (Peebles (1980)). 

We can note th at this result agree s with the calcula- 
tions performed in Bharadwaj (1996) from the BBGKY 
hierarchy. Indeed, in that paper it was shown that the 
lowest order non-linear correction to the two-point corre- 
lation function derived from the BBGKY hierarchy (which 
remains valid in the non-linear regime after shell-crossing) 
matches the prediction of the standard hydrodynamical 
approach. Note that the BBGKY hierarchy can be de- 
rived from the coUisionless Boltzmann equation (yj), see 
Peebles (1980) . Therefore, both works explicitly show that 
multi-streaming effects cannot be handled through a per- 
turbative method. 



4.2.2. Peculiar velocity 

From the distribution function Sf we can also derive the 
mean velocity of the fiuid. Indeed, the hydrodynamics 
equations are obtained from moments of the Boltzmann 
equation. In particular, the hydrodynamic peculiar veloc- 
ity is defined by: 



pv= dp/(x,p,T) -. 
Since p = p{l + S) we obtain using eq.(Eh: 
[1 + (5(x, r)] v(x, r) = e^^r J ^p Sf{^, p, r) p 
which yields for the second-order term: 



(49) 
(50) 

(51) 

"2'^Pi which 

is given by eq. (E9) . Then, we obtain the second-order term 
from eq.(E3). This yields: 



.(2) 



dp Sf 



(2) 



5(i)v(i 



The first-order term is simply v 



(1) 



VL 



,(2) 



(k, 



-4ki 
I5 fcf 



" /"dki,5Lo(ki)5Lo(k-ki) 
4 k.(k-ki)ki 6 k.ki k 



15 Ik 



ki|2 fc2 35 kl P 



4 k.kik.(k-ki) k 

'35~fcf |k-ki|2 fc2 



(52) 
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Within the hydrodynamical approach, when one expands 
the density and velocity fields over the linear growing 
mode (i.e. there are no decaying modes) the velocity field 
is nonrotational (e.g., Peebles (1980)| ). Then it is fully de- 
fined by its divergence and one usually introduces: 

6l(x, r) = ^, hence e{k, r) = - e^ ik.v(k, r). (53) 

a 2 

In real space, using eq.( p2|) we recover th e result of the 
hydrodynamical approach ( Peebles (1980) ): 



0(2)(x,t) =-e4^ydkidk2 e'(''^+''^)-^<5Lo(ki)<5Lo(k2) 



3 ki.ka 4 /ki.ka 
7 kj 7 V kik2 



(54) 



In a similar fashion, we can consider the normalized vor- 
ticity u; = (V X v)/d. Using eq.(p2) we can easily check 
that we obtain uj'^^^ — (and of course the linear mode as 
defined in eq.(En) is nonrotational). Thus, we see that up 
to the second-order we recover the results of the hydrody- 
namical approach. 

5. Divergence of perturbative series 

In the previous section we noticed that the perturbative 
treatment of the coUisionless Boltzmann equation yields 
the same results at second-order as the hydrodynamical 
approach when the linear mode is given by eq.(E3). In 
fact, it is easy to see that this equality extends to all 
orders of perturbation theory. Indeed, let us consider an 



initial condition of the form (43) which is very smooth. 
That is, (5Lo(k) shows no power at small scales (i.e. large 
wavenumbers k) and the flow is nonrotational, defined by 
the linear growing mode of hydrodynamics. Then, there 
is no shell-crossing until a finite time Tc (i.e. tc > 0). 
This means that the hydrodynamical description is ac- 
tually exact at early times t < tc'- at each point there 
is only one velocity. Then, until the non-zero time tc the 
solution of hydrodynamics (i.e. continuity and Euler equa- 
tions) is also a solution of the coUisionless Boltzmann 
equation from which they derive. This implies that the 
density contrast and velocity fields obtained by the pertur- 
bative theory in both frameworks are identical. Moreover, 
the form of the perturbative expansion does not depend 
on the properties of (^Lo(k), hence this identity holds for 
any initial conditions Slo{^)- As a consequence, if the lin- 
ear mode is of the form of eq.(^ the perturbative results 
obtained from the coUisionless Boltzmann equation and 
from the hydrodynamical description are identical. Here, 
we note that this identity was explicitly derived at all 
orders within the framework of the Zel'dovich approxi- 
mation in Bharadwaj (1996)| . Therefore, in order to study 
the non-linear regime one needs to devise non-perturbative 
methods applied to the coUisionless Boltzmann equation 

(0). 

On the other hand, in the case of CDM-like power- 
spectra (or power-law power-spectra) there is some power 



at all scales. In particular, the density fluctuations increase 
with no limit at small-scales (if we do not include a cutoff 
at small scales). Then, as soon as i > there is some shell- 
crossing. This implies that the hydrodynamical solution is 
no longer a solution of the coUisionless Boltzmann equa- 
tion (actually, the former is not well deflned). Then, the 
perturbative expansion must diverge (otherwise the hy- 
drodynamical solution and the Boltzmann solution would 
be identical since they would be equal to the same per- 
turbative solution). Hence, we conclude that for hierar- 
chical initial conditions the perturbative expansion is only 
asymptotic. If we include a cutoff at small scales for the 
initial power-spectrum the perturbative series may con- 
verge until the finite time tc when shell-crossing occurs 
for the first time (though this is not guaranteed a priori). 
However, it must diverge at later times which implies that 
the perturbative series diverges at all times of practical in- 
terest. 

In fact, the system we investigate here can actually be 
even more pathological than these results suggest. First, 
as we show in a companion paper (paper II) where we 
develop a non-perturbative method to study the quasi- 
linear regime, the generating function of the cumulants of 
the density contrast obtained at leading order has a ra- 
dius of convergence which is zero independently of shell- 
crossing for P{k) ex fc" with n < 0. This means that the 
high-density tail of the probability distribution function 
of the density contrast cannot be evaluated by perturba- 
tive means in this case. Secondly (and more importantly), 
the terms encountered in perturbative expansions actually 
diverge beyond some finite order for power-law power- 
spectra (e.g., Scoccimarro & Frieman (1996), paper V). 



This implies that one can only use perturbation theory 
up to a finite order at best. 

6. Other cosmologies 

In this section, we show how we can extend the perturba- 
tive calculations developed in Sect.|| for a critical density 
universe to arbitrary cosmological parameters. In the gen- 
eral case, the scale factor a{t) is no longer a power-law 
hence it is not useful to introduce the time coordinate r 
as in eq. (pTJ) . Thus, from eq.(^ and keeping the coordinate 
t we define the operator O by: 



(0.,5/)(k,p,t) 



dSf k.p 

ot a-' 



AnQp k dSp 
a fc2 dp 



dp'5f{k,p',t). 



(55) 



Next, as in Sect . [4. l| we obtain the growing mode r]{'k, p, t) 
by keeping only the linear terms of the linear hydrody- 
namical growing mode which can still be written as in 
eq.M). This yields: 



V = D+SLoi^)SDip) 



.a^I)+5,o(k)^.^(p) (56) 



where D^{t) is the usual linear growing mode, which is 
no longer proportional to a{t). Then, one can easily check 
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that O.rj = with 77 given by e g. ([5^) and O by eq.(| 
using the fact that D+{t) obeys (iPceijles (1980)1): 



D^ 



2^D, 



inGp 



D. 



(57) 



Then, in order to get the second-order term we simply keep 
the structure of the result (El) obtained for the critical 
density universe but we replace all time-dependent factors 
by unknown functions. Thus, we write: 

^/(2) (k, P,t)^ J dki ,5Lo(ki)5Lo(k - ki) 
f ^•'^ix r \ , k.kik.(k-ki) 



kf 



K |k-ki| 



. ki dSo , . k.(k-ki)ki dSo , . 

"'^^fc? -^^P) " 'g^ |k-kiP fcf-^^P^ 

. k.ki k dSo, , k.kik.(k-ki) \l dSo, , 

-*55^^r2 ■7^(P) - »g6 ,2 „. ,. I. — -^ztIp) 



kl fc2 ■ ap 



/c2 |k-ki|2 p- ap 



k]_ _5_ / k~ki 9(5g 
"^'fc^ap ^|k-ki|2- ap ^P^ 



(58) 



where gi{t), ...,g-j{t) are functions of time. Next, we sub- 
stitute eq. (^ into eq. (^ , using eq. (|6|) for the first-order 
term of (5/, and keeping only the second-order terms we 
obtain a system of differential equations for the unknown 
functions gi{t), ...,g7{t). This yields: 

a^ffi = 53 + 55 , a^52 = 54 + 56 , 053 = AirOpDl (59) 
a^54 = 2^7 , 055 = 47r^p5ri (60) 

age = 4,TrGpg2 , 57 = 47rC/pi:)^ (61) 

which determines the functions 51 (t), ...,g7{t). The bound- 
ary condition at i = is given by the constraint that for 
i ^ we must recover the result (^) obtained for the crit- 
ical density universe. We shall not go further here since 
the discussion of Sect.g remains valid for arbitrary cos- 
mologies hence we already know that we shall eventually 
recover the results obtained through the hydrodynamical 



approach, described in Bouchet et al. (1992) 



Thus, in order to derive the perturbative expansion of 
the distribution function df in the general case one first 
computes the terms Sf^"' up to the required order for 
a critical density universe, using the method developed 
in Sect|3[ This provides an explicit expression for 5/^"^ 
Then, one looks for a solution of eq.(Q) in the case of in- 
terest (which differs from the critical density universe) by 
keeping the same structure for (S/^") but replacing all time- 
dependent factors by unknown functions of time. Then, 
these functions are determined by substitution into eq.(0). 
The need for this two-step procedure comes from the fact 
that in the general case the calculation of the operator 
C_i is not straightforward. More precisely, in order to 
derive O-i one can st ill follow the procedure outlined in 
App.^ up to eq.( A.12| ). However, the recursion obeyed by 
the functions i?,i(cr) is no longer a simple Laplace convo- 
lution so that it is not obvious to find a simple analytic 
solution. 



7. Functional formulation 

Finally, we briefiy describe how one can use the integral 
equation ( |3l| ) to obtain an explicit expression for the sta- 
tistical properties of the distribution function. In this sec- 
tion, we shall work in real space x, and we note u: the 
7-dimensional coordinate u: — (x, p, r). Thus, the fields 
77(0;) and Sf{ijj) are real. Taking the Fourier transform 
of eq.(pl|) we see that 5f{u}) and 77(0;) are related by the 
same quadratic equation where the kernel iir(r;ri,r2) is 
replaced by its real-space Fourier transform K{u};u}i,uj2)- 
In practice we are not really interested in the exact solu- 
tion 6f associated with a given field rj. Indeed, since the 
initial condition 77 is a random field we are only interested 
in the statistical properties of the distribution function 6f. 
These can be derived from the functional Z[j] of the test 
field j{u}) defined by: 

Z[j]^{J'"^^-'f) (62) 

where (..) expresses the average over the initial conditions. 
Expanding the exponential we get: 

Z[j] = 1 + X!~i / du;i..da;„ j(a;i)..j(u;„) 

ri=l ^' •' 

X (,5/(a;i)..,5/(a;„)) (63) 

which clearly shows that Z[j] yields all moments 
{df{uJi)..Sf{uJn)) of the distribution function Sf. If we 
assume the initial conditions to be Gaussian we can write 
the average (p2[) as the path-integral: 



Z[j]^Mo / [d77(a;)] e^-*^'-HA-^'7 



(64) 



where we introduced the normalization constant: 

Afo^{BetA-'Y^\ (65) 

Here 5f[ri] must be understood as the solution of eq.(^ 
(in real w-space) associated with a given 77. We also used 
the short-hand notations: j.Sf = J du) j{u)).Sf{u)) and 
77.A,Y^?7 = /da;ida;2 77(a;i).A,yi(a;i,a;2). 77(0^2). The ker- 
nel A^^ is the inverse of the kernel A^ defined by: 

A^(a;i,a;2) = (77(u;i)77(u;2)) (66) 

which fully determines the statistics of the random field 
77(0;) since the latter is assumed to be Gaussian. Here we 
must note that the kernel A^ is not invertible. This is 
related to the fact that the linear mode 77(0;) is restricted 
to the form (MS), which also ensures that O.rj — 0. Thus, 
in eq.(p4) it is understood that the kernels A,, and A~^ 
are regularized through a small modification described by 
an infinitesimal parameter e, which makes A,, and A""'^ 
invertible. For instance, one can add to the kernel A^ a 
term of the form e x (5£)(u;i — 0^2). Then, the final results 
are obtained by taking the limit e — > in the calculations. 
Using eq. ( |3l| ) we can make the change of variable 77 —> 5f 
in the expression (^4|). This yields: 

Z[j] = AAoy[d,5/(a;)]|DetM| 
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where we used the notation: 

{KS f)iu:)= I du:^<:\u^2K{u:-u:i,uJ2)5f{uJi)5f{uJ2)- (68) 

In eq.(|67|) we introduced the matrix M defined by: 

5{df[u:2)) 
-<5B(a;i-u;2)-2.g j dJ k,{u,i;u>2,u,')6f{u,'). (69) 



Here we used eq. ( |3l| ) (in real w-space) and Kg is the sym- 
metrized part of K, as in eq.(p2|). Next, we need to com- 
pute the Jacobian |DetM|. To do so, we use the relation: 



DctA/^e^-^'"^^ 

and from eq. (|6S|) we have: 

Tr In A/ 



VM:TrA" 



(70) 



(71) 



n=l 



where we introduced the matrix A: 

A(a;i, u;2) = j du;' k,{uj^-u:2, uj') df{io'). (72) 

Hence a simple recursion yields for n > 2: 



Next, we note from eq.(B.l) that the kernel Ks{u}; u}i,u^2) 
involves a factor 6{t — Ti)(5£)(r2 — ti) (which also ensures 
that causality is satisfied) . Then, we see that the integrand 
in eq.([73|) involves a factor 6'(r— ri)..0(T„_i— r'). Thus, the 
Heaviside factors imply that the time-coordinates must 
obey: r > ti > .. > r„_i > r' in order that the integrand 
be non-zero. This also implies that v4"(a;,a;) = for n > 
2. On the other hand, the trace TrA" is defined by: 



TrA"= /da;A"(a;,a;). 



(74) 



Hence we see that TrA" = for n > 2. Note that 
this result only relies on causality requirements (which 
translate into the Heaviside factors) and it also applies 
to usual stoc hastic differential e quations (e.g., Langevin 
equation, see Zinn-Justin (1989) ). Thus, we are left with 
DetAf = exp(-2gTr^), with: 

2TrA := 2 / dwdw' Ks{uj;lj,lj') Sf{uj') 



K{u}-,u},u;') + K{u;;lj',u}) Sf{u;'). (75) 



= / du;da;' 



Then, from cq.( P3.l|) we can see that: 

du; K{u:; u:, a;') - / da; X(a;; a;', a;) = 0. 



Indeed, the second term in eq.(BJ) vanishes since t — ti 
and we can push the integration path over si to the right: 



Re(si) — > oo. The first term vanishes after we integrate 
over si and next over p or over the angles of k' (using the 
fact that K{u}; w', uj) oc (5£i(k')). Note that this latter inte- 
gration involves the "regularization" of the gravitational 
interaction discussed in Sect.0 below eq.(P). Indeed, as 
discussed in Sectfl apparently undetermined integrals are 
given a meaning by the prescription that one must first 
perform the integration over angles (or introduce a large- 
scale cutoff kc — > 0). Thus, we get DetM = 1 and we can 
write eq.(|67|) as: 

Z[j]=Mo /'[d<5/(u;)]e^-^/-^^/-^.-^-^/+*^35/^-4if45/- 

(77) 
Here we used the short-hand notations: 

K3 Sf = XI / da;da;i..da;3 6f{uJi)6f{uj2)Sf{uj3) 
perm. "^ 

X A-\u:i,u:)K{u;u;2,u:3) (78) 

and: 

KiSf^ = 2 X! / du;du;'du;i..du;4 Sf{u}i)Sf{u}2)Sf{u}3) 
perm. -^ 

X Sfiu;^)kiu;;u:,,u:2)^^H^,^')Kiu:';uj3,^4) (79) 

where the sums run over the permutations of (wi, a;2, ^3) 
and (wi, ci;2, ^3, ^4). We symmetrized the kernels K3 and 
ii'4 in eq. (fTq) and eq.(u9) because it simplifies perturba- 
tive expansions. Thus, the path-integral ( |77| ) yields an 
explicit expression for the functional Z[j], hence for all 
statistical properties of the distribution function 5f{uj). 
Indeed, all terms in this expression are kno wn. The kernel 
K which yields A'3 and K^ is given in eq.(BT) while the 



kernel A is obtained from A^ (after a suitable regular- 



ization). The latter is given by eq. (|66|) and eq. ([43|) and it 
can be expressed in terms of the power-spectrum P{k) of 
the linear density field 5n){k). 

The correlation functions can be obtained from eq. ( [77| ) 
as a perturbative series over the coupling constant g. 
Indeed, we noticed in Sect. |3.3| that this is equivalent to the 
usual expansion over the linear growing mode ?y(u>). To do 
so we simply develop the exponential term in eq.(77) over 
gK-i and g^K^ as a series over g. This corresponds exactly 
to the Feynman diagrams in standard Quantum Field 
Theories for a "(/-^ + 0''" theory (e.g., ^inn- Justin (1989)D 
(note that we defined the kernels K3 and /Q in such a way 
that they are symmetric) . One must merely pay attention 
to the fact that the kernels K3 and K^ are not pure num- 
bers. Of course, this is due to the long-range character 
of gravity which leads to a non-local theory (while usual 
Quantum Field Theories are local). Note that we do not 
need the explicit expression of the inverse kernel A~^. For 
instance, in order to compute the skewness some terms 
(76) which depend on A~^ appear but they actually cancel 
one another. 

In practice, the expression ( [77| ) is not very conve- 
nient for perturbative calculations and it is easier to use 
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the direct method described in Sect.0, or simply use the 
standard hydrodynamical approach. However, we stress 
that eq. I jn\) presents the advantage to be an explicit non- 
perturbative expression for the functional Z[j]. Thus, one 
may hope to be able to extract some useful information 
from this relation for the non-linear regime. However, this 
is beyond the scope of this article. 

8. Conclusion 



Thus, in this article we have developed a system- 
atic method to obtain the solution of the collisionless 
Boltzmann equation, which describes large-scale struc- 
ture formation in the universe, as a perturbative series. 
Moreover, we have explained that these perturbative re- 
sults are identical to those obtained from the hydrody- 
namical description of the fluid, if we start with the same 
initial conditions. Thus, multi-streaming cannot be han- 
dled through perturbative methods. Besides, for hierar- 
chical scenarios the perturbative expansions diverge. The 
interest of this approach is that, contrary to the hydrody- 
namical model, the collisionless Boltzmann equation pro- 
vides a rigorous description of the dynamics, even in the 
non-linear regime. Hence in order to tackle the non-linear 
regime where shell-crossing plays a key role one should 
use these equations of motions. Then, the results pre- 
sented in this paper may serve as a basis for a study of the 
non-linear regime. In particular, we have obtained a non- 
perturbative integral equation (|3l|) which explicitly relates 
the non-linear distribution function to the initial condi- 
tions. Besides, we have derived an explicit path-integral 
expression (^^ which yields the correlation functions of 
the actual non- linear density field at all orders. Although 
this result provides a formal description of the statis- 
tics of the density field one still needs to check whether 
it is possible to extract some practical information from 
this expression. Thus, we shall describe in companion pa- 
pers an alternative non-perturbative method, based on the 
path- integral (K4), which allows one to derive rigorous re- 
sults in the quasi-linear regime for the probability distri- 
bution of the density contrast for Gaussian initial condi- 
tions (paper II) as well as for some non-Gaussian scenarios 
(paper III). Such a formalism can also be used in the non- 
linear regime (paper IV). 



Appendix A: Calculation of the exponential e "^^ 

In order to obtain the kernel K defined in cq. (psj ) we first 
need to compute the exponential e^°''-' . Thus, for a given 
function F{r) = i^(k, p, r) we look for the function B{r, a) 
defined by: 



B(r,fT)EEe-^'^M.i^ 



(A.l) 



To derive B{r, a) we note that this function satisfies the 
initial-value problem: 



dB 



-O.B and B(r,0)=F(r). 



Moreover, we are only interested in the restriction of 
B{v,a) to (7 > 0. Using eq.(Qj) we obtain the first-order 
integro-differential equation: 

dB dB , , ^ 

-3i(k.p) e-^ B 



da Ot 



-2ie 



r k dSp 
dp 



{p)Jdp'B{k,p\r,<j). (A.3) 

This equation is formally solved by the method of charac- 
teristics which yields the integral form: 

S(k, p, T + a, a) = B(k, p, T, oy'^ii^-p)---[e--^] 



k dSr 



-2z^.^^(p) / da'e^^^'^P)^""! 
/o 



fc2 dp 
X e^^"' I dp' B(k, p', r + <t', a') 



(A.4) 



The form of eq.(A.4) shows that we can consider in the 
next intermediate steps k and r to be fixed constants and 
we are led to introduce the functions: 



G(p, a) = B{k, p, r + cr, a) 
and: 

C/(P,(t) = B(k,p,T,0) e3'(k.p)e-[e---l] 

as well as the kernel: 

JC{p,a\p\a') ^ 9{a-<j')^.^{p)e'^ 



(A.5) 
(A.6) 



X e 



3i(k.p)e"^[e"''-e 



(A.7) 



where 9{a — a') is Heaviside's function. Then, eq.( |A.4| ) 
reads: 

G(p,a) = Uip,a)+iy / /C(p, a|p', a')G(p', (T')dp'da'(A.8) 

where we introduced the constant: 

i^ = -2i e^ (A.9) 

and we integrate over a' from to -l-oo: J da' = L da' . 
The linear integral equation ( A. 8 ) is a standard Fredholm 
equation of the second kind and it can be solved as a 
perturbative series over v. 

G{p,a)^U{p,a) 

+ Zl ''" / ^"(P' '^Ip'' ^')t^(p', cr')dp'da' (A.IO) 

n=l -^ 

where the kernels /C„ are defined by the recursion: 

/C„ = y dp"da"/C(p,a|p",a")/C„_i(p",a"|p',(7') (A.ll) 

with /Ci(p,(t|p', cr') = IC{p,a\p',a'). Thus, in order to 
obtain G(p, a) we simply need to sum up the scries in 
the r.h.s. of eq.(AT^). One can check by substitution into 
eq.(A.ll|) that these kernels /C„ are given for n > 2 by: 



k dSo , w „. -r\"-l 



(A.2) 



X e 



3i(k.p) 



dcr„_i 9{a - an-i) 

'1 Rn{cTn-l) (A.12) 
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where we introduced: 



R7i{o'n-l) — 



X e 



dcri..dCT„_2 ^(o-„-i - cr„-2) 

,-1 [e-'^-i - e"-'] 



(CTi -o- 



function tpn (here we add a subscript "K" to distinguish 
Kummer's function ipx from ip defined in eq.(A.17)) by 
(e.g., Batenian (1953)): 



■0(5, x) — e^x*r[— s, x] — "0^(1, 1 — s; x). 



(A.19) 



for 71 > 3 and i?2(ci) = 0{ai—a') 



(A. 13) This implies that 'ip{s,x) is an entire function of s. 

In particular, one obtains from eq.(A.19) the asymp- 

One can see totic behaviour for large negative s at fixed x (see 



from eg. ( A. 13 ) that the functions Rn{cr) actually satisfy 
the recursion: 



Bateman (1953)| , vol.I,§6.13.2.(6)): 

-oo : iIj{s, x) 



Re(s) 



2,SgS-(s+l/2)ln(-s) 



(A.20) 



Rn{<j) 



dl 1 



„<y-i\ 



Rn-l{l) 



(A. 14) which shows that ijjjs^x) diver ges for Re(s) — > — oo. On 
the other hand, from eq.( A.17 ) we see that the function 
which has the form of a simple Laplace convolution. This i^{s,x) also satisfies the relation (inverse Laplace trans- 
recursion is easily solved by taking the Laplace trans- form): 
form of eg. (A. 14) which yields, after going back to cr-space 



through an inverse Laplace transform: 



-siij+x[l — e'^ 



Rn{a) 



+ 100 



ds 
27ri 



s{cr — cr') 



-1 



-\-ioo 



ds_ 
27ri 



e"" il;{s + si,x). 



(A.21) 



(A.ioj Then, we can write eq.( A.16 ) as: 



Here and in the following calculations the integration path 
over the Laplace variables like s will always be to the right 
of the singularities of the integrand. Thus, in eg.( A.K) the 
integration path for s obeys: Re(s) > 1. Then, eg. ( A. 15 ) 
provides an explicit expression for the kernels /C„, see 
eq.(A.12). Next, we can now compute the sum over n in 
eq.(A.lO); which eventually provides an explicit expres- 
sion for B{r, a) using eg.( A^). Thus, a simple calculation 
gives: 



S(k,p,T,a)= / ^e^^^^(si,3z(k.p)e--)^F(k,p,r-a) 



^ k dSo , . 



ds2ds3 



,(s2 + S3)o-. 



Si(si - 1) 



(2^1)2-- (si-3)(.si+2) 

X —-^ [ dp' ^{s2,-3i{k.p')e--) 

Sl + S2-1 J 

X ^{S3 + l,3z(k.p')e-") F(k,p',T - a)}. (A.22) 



B{k, p, T, a) = e 
k dSo 



-2i- 



ds ,./„ „i 



fc2 dp 

„ „T+cr'-a ^{^ 



ip) J dp' da' da, J 



Using the definition (^2|) this result provides the operator 
O-i through the formal integral: 



2ni 



cr-CTi 



]+3J(k.p')e""[e 



(0_i.F)(r)= / daB{r,a). 



1) 



(s-3)(s-l-2) 



F{k,p',T-a). 



(A.16) 



(A.23) 



Appendix B: Calculation of the kernel ^(r, ri,r2) 



It is conv enien t to replace the integrations over a' and From eg.( |A.22D and the expression ( |A.23| ) for the oper- 
(Ji in eg.( |A.16| ) by integrals over the conjugate Laplace ator 0_i we can derive the explicit form of the kernel 
coordinates s' and si. To do so, we simply compute the K{r;ri,V2) defined in eg.(|25|), using eg.(|l|). First, we 



ipis,x) 



da e 



— sa-\-x[l — e^ 



integrals over a' and cti in eg.( |A.16D in terms of the func- write the factor SoiTi - t)(5_d(t2 - r) which a ppear s in 

tion ^p{s, x) which we define by: eg.(|l|) as 6d{ti - t)5d{t2 - n). Then, using eg.( |A.22| ) we 

see that -B(r, a) — e~'^*-' .K exhibits a factor SoiTi — T+a). 

(A 17) Thus, the integration on a over the range [0,oo[ is now 

/o straightforward and it yields a Heaviside factor 6{t — ti). 

In the following the argument x of the function ^{s, x) will ^ext, a simple calculation yields for k{v- n, r2): 

be imaginary. Then, we have the asymptotic behaviour for k ^2i e"^ 0(t - ti) SdW^ ~ n) ^^(ki -I- k2 - k) 
large positive s: „ , r i nx 



Re(a;) = 0, Re(s) -^ -|-oo : i^isjx) 



(A.18) 



since for imaginary x and Re(s) > eq.(A.17) i mplies : 
1-0(5, a;)| < 1/Re(s). The integral representation ( A. 17 ) 
only applies to Re(s) > however the fun ction can be 
continued to Re(s) < 0. Indeed, from eg. ( A. 17 ) one can 
check that ^{s^x) can be expressed in terms of the incom- 
plete Gamma function F or the confluent hypergeometric 



ds2ds 



3 g(s2+S3)(T-ri) 



ki.k k_ d5D_ 

kf fc2 ■ 5p (P) J (27ri)2 
(Sl - 3)(S1 + 2)(S1 + S2)(S1 + S2 - 1) 

X V(s2, -3i(k.p2)e"^)-(/'(s3, 3i(k.p2)e"^; 



(B.l) 
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where the integration paths for the variables si and S2 
obey: Re(si) > 3 and Re(s2) > 0. To get eq.(B.l) we used 
the following properties oi ^(s, x) which are obvious from 
eq.(|Al7|): 



^(s,0) 



dx 



(s, x) = ip{s, x) — ip{s — 1, x). 



(B.2) 



Note the factor 0{t — ti) in eq.(BJ_) which explicitly shows 
that the expansion ( |2l| ) satisfies caus ality requirements. 



This agrees with the discussion in Sect. 3.2 



Finally, as discussed in Sect. 3.2 we need to check 
that the kernel ir(r;ri,r2) obtained in eq.(B.l) satisfies 
eg. (20 ). To do so, it is convenient to first integrate over si 
in eq.(B.l). The integration of the first term is straight- 
forward, using the inverse Laplace transform (A.21). In 
order to integrate the second term over si we use the iden- 
tity (pq). This allows us to express this factor in terms of 
ip{si,0) — 1/si and iIj{si — 1,0) = l/(si — 1), so that the 
integration over si becomes quite simple. In particular, 
the contribution from the poles si = —S2 and si = 1 — S2 
vanishes when we take the integration path over S2 to be 
Re(s2) > 3. Then, we can write the kernel K as: 

K^2t e"i Sd{t2 - Ti) 5i5(ki + k2 - k) 



6 ki.k 



+«oo 



dS2dS3 



,(s2+S3)(t-Ti) 



5 kl J_,^ (2^j)2 

X V(s2,-3i(k.p2)e"'') -0 (s3,3i(k.p2)e"^) 



k dS 



D 



-(P) 

+3ie"^(5_D(p) 



2 e3(T-ri) 



3 e-2(— -i) 



fc2 • dp '^' V (S2 + 2)(S2 + 3) (S2 - 2)(S2 - 3) 

eHr-ri) p-2(r-ri) 



(s2 + 2)(s2 + 3) (s2-2)(s2-3) 

(B.3) 

where the integration path obeys Re(s2) > 3. Note that 
we removed the Heaviside factor 9(t — ti) from the second 
term in eq.(B.3) because the integrals over S2 and S3 van- 
ish for (r — Ti ) < (since for (t — ri ) < we can push the 
integration path to Re(s) — > -|-cxi). Next, we can apply the 
operator O defined in eq.([l4) onto the kernel K written 
in eq.(B.3). After a tedious calculation we finally obtain 
eq.(20). The only technical trick in this calculation is to 
use eq.(B.2) as well as the property: 



s ip{s,x) = 1 — X tpis — 1, a;). 



(B.4) 



This identity is obtained by multiplying eq.( |A.17 ) by s 
and integrating by parts. This explicit check of eq.(p0[) 
fully justifies the procedure developed in Sectjsj 
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